#ifndef LB_Sparse_Matrix_H_
#define LB_Sparse_Matrix_H_
#define Sparse_Matrix LB_Sparse_Matrix 
#include<Math/LB_Math.h>

#ifdef __cplusplus
extern "C" {
#endif

// typedef struct LB_Sparse_Matrixn{
//     int rows,cols;
//     int*i_index;
//     int n;
//     double**data;
// }LB_Sparse_Matrixn;
// static inline void lb_sparse_matrixn_init(LB_Sparse_Matrixn*lsm)
// {
//     lsm->i_index=NULL;
//     lsm->data=NULL;
//     lsm->rows=0;lsm->cols=0;lsm->n=0;
// }
// static inline void lb_sparse_matrixn_reserve(LB_Sparse_Matrixn*lsm,int n)
// {
//     lsm->i_index=(int*)malloc(sizeof(int)*n);
//     memset(lsm->i_index,0,sizeof(int)*n);

//     lsm->n=n;
//     lsm->data=(double**)malloc(sizeof(double*)*lsm->rows);
//     for(int i=0;i<n;i++)
//     {
//         lsm->data[i]=(double*)malloc(sizeof(double)*n);
//         memset(lsm->data[i],0,sizeof(double)*n);
//     } 
// }










typedef struct LB_Sparse_Matrix_Node{
    int i,j;
    double data;
    struct LB_Sparse_Matrix_Node* right,*down;
    void* prop;
}LB_Sparse_Matrix_Node;
static inline void lb_sparse_matrix_node_init(LB_Sparse_Matrix_Node* lsmn )
{   
    lsmn->i=0;lsmn->j=0;lsmn->data=0;
    lsmn->right=NULL;lsmn->down=NULL;
    lsmn->prop=NULL;
}
static inline void lb_sparse_matrix_node_free(LB_Sparse_Matrix_Node* lsmn )
{   
    SAFE_FREE(lsmn);
}
typedef struct LB_Sparse_Matrix{
    LB_Sparse_Matrix_Node** rhead,**chead;
    int rows,cols,n_nz;
    void *prop;
}LB_Sparse_Matrix;

static inline void lb_sparse_matrix_init(LB_Sparse_Matrix*lsm)
{
    lsm->rhead=NULL;lsm->chead=NULL;
    lsm->rows=0;lsm->cols=0;lsm->n_nz=0;
    lsm->prop=NULL;
}
void lb_sparse_matrix_free(LB_Sparse_Matrix*lsm);



static inline void lb_sparse_matrix_initn(LB_Sparse_Matrix*lsm,int rows,int cols )
{
    lsm->rows=rows;lsm->cols=cols;lsm->n_nz=0;
    lsm->prop=NULL;
    lsm->rhead=(LB_Sparse_Matrix_Node**)malloc(sizeof(LB_Sparse_Matrix_Node*)*(lsm->rows+1));
    lsm->chead=(LB_Sparse_Matrix_Node**)malloc(sizeof(LB_Sparse_Matrix_Node*)*(lsm->cols+1));
    
    memset(lsm->rhead,0,sizeof(LB_Sparse_Matrix_Node*)*(lsm->rows+1));
    memset(lsm->chead,0,sizeof(LB_Sparse_Matrix_Node*)*(lsm->cols+1));
}

void lb_sparse_matrix_insert(LB_Sparse_Matrix*lsm,int i,int j,double data);




static inline double lb_sparse_matrix_get_data(LB_Sparse_Matrix*lsm,int i,int j)
{
    LB_Sparse_Matrix_Node*q=NULL; 
    for(q=lsm->rhead[i];q!=NULL&&q->j<j;q=q->right);
    if(q==NULL||q->j!=j)
    {
        return 0;
    }
    return q->data;
}


// static inline double lb_sparse_matrix_ij_mult(LB_Sparse_Matrix_Node* rhead,LB_Sparse_Matrix_Node* chead)
// {
//     for(int i=0;)
//     {



//     }

// }

static inline LB_Sparse_Matrix* lb_sparse_matrix_mult_sparse_matrix(LB_Sparse_Matrix*lsm1,LB_Sparse_Matrix*lsm2)
{
    if(lsm1->cols!=lsm2->rows){return NULL;}
    LB_Sparse_Matrix* re=(LB_Sparse_Matrix*)malloc(sizeof(LB_Sparse_Matrix));
    lb_sparse_matrix_initn(re,lsm1->rows,lsm2->cols);
    //lb_sparse_matrix_init(re);
    LB_Sparse_Matrix_Node*q=NULL,*p=NULL;  
    for(int i=0;i<lsm1->rows;i++)
    {
        for(int j=0;j<lsm2->cols;j++)
        {
            double temp=0;
            p=lsm2->chead[j];
            for(q=lsm1->rhead[i];q!=NULL;q=q->right)
            {

                for(;p!=NULL&&p->i<q->j;p=p->down);
                if(p==NULL)
                {break;}
                if(p->i==q->j)
                {
                    temp+=p->data*q->data;
                }  
            }
            if(fabs(temp)>1e-9)
            {
                lb_sparse_matrix_insert(re,i,j,temp);
            }
        }
    }
    return re;
}

static inline double* lb_sparse_matrix_mult_array(LB_Sparse_Matrix* lsm,double* vec)
{
    double *re=(double*)malloc(sizeof(double)*lsm->rows);
    memset(re,0,sizeof(double)*lsm->rows);
    LB_Sparse_Matrix_Node*q=NULL;
    for(int i=0;i<lsm->rows;i++)
    {
        for(q=lsm->rhead[i];q!=NULL;q=q->right)
        {
            re[i]+=q->data*vec[q->j];
        }
    }
    return re;
}


// 
// 要保证lsm是对陈正定矩阵
//
double* compute_conjugate_gradient(LB_Sparse_Matrix*lsm,double* b);

// 
// 利用prop属性
//
static inline LB_Sparse_Matrix* lb_sparse_matrix_transepose(LB_Sparse_Matrix*lsm)
{
    //printf("here %d\n",lsm->n_nz);
    if(lsm==NULL||lsm->n_nz<=0){return NULL;}
    void** props=(void**)malloc(sizeof(void*)*lsm->n_nz);
    int num=0;
    LB_Sparse_Matrix_Node*q=NULL;
    for(int i=0;i<lsm->rows;i++)
    {
        for(q=lsm->rhead[i];q!=NULL;q=q->right)
        {
            // re[i]+=q->data*vec[q->j];
            props[num]=q->prop;
            LB_Sparse_Matrix_Node* p=
                (LB_Sparse_Matrix_Node*)malloc(sizeof(LB_Sparse_Matrix_Node));
            lb_sparse_matrix_node_init(p);    
            p->i=q->j;p->j=q->i;p->data=q->data;
            q->prop=p;
            num++;
        }
    }
    LB_Sparse_Matrix* re=(LB_Sparse_Matrix*)malloc(sizeof(LB_Sparse_Matrix));
    lb_sparse_matrix_initn(re,lsm->cols,lsm->rows);
    re->n_nz=lsm->n_nz;
    for(int i=0;i<lsm->rows;i++)
    {
        q=lsm->rhead[i];
        LB_Sparse_Matrix_Node* p=(LB_Sparse_Matrix_Node*)(q->prop);
        re->chead[i]=p;
        for(;q!=NULL;q=q->right)
        {
            p=(LB_Sparse_Matrix_Node*)(q->prop);
            if(q->right!=NULL)
            {
                p->down=(LB_Sparse_Matrix_Node*)(q->right->prop);
            }
        }
    }
    
    for(int i=0;i<lsm->cols;i++)
    {
        q=lsm->chead[i];
        LB_Sparse_Matrix_Node* p=(LB_Sparse_Matrix_Node*)(q->prop);
        re->rhead[i]=p;
        for(;q!=NULL;q=q->down)
        {
            p=(LB_Sparse_Matrix_Node*)(q->prop);
            if(q->down!=NULL)
            {
                p->right=(LB_Sparse_Matrix_Node*)(q->down->prop);
            }
        }
    }

    //printf("here %d\n",lsm->n_nz);
    num=0;
    for(int i=0;i<lsm->rows;i++)
    {
        for(q=lsm->rhead[i];q!=NULL;q=q->right)
        {
            // re[i]+=q->data*vec[q->j];
            q->prop=props[num];
            num++;
        }
    }
    free(props);
    return re;
}


static inline void lb_sparse_matrix_printself(LB_Sparse_Matrix*lsm)
{
     if (lsm == NULL) {
        printf("Sparse Matrix is NULL.\n");
        return;
    }

    printf("Sparse Matrix:\n");
    printf("Rows: %d, Cols: %d, Non-Zero Elements: %d\n", lsm->rows, lsm->cols, lsm->n_nz);

    // 遍历每一行
    for (int i = 0; i < lsm->rows; i++) {
        LB_Sparse_Matrix_Node *current_row = lsm->rhead[i];
        if(current_row==NULL)
        {
            continue;
        }
        // 在当前行中遍历每一个非零元素
        while (current_row != NULL) {
            printf("(%d, %d ,%lf) ", current_row->i, current_row->j, current_row->data);
            current_row = current_row->right;
        }
        printf("\n");
    }
}


#ifdef __cplusplus
}
#endif

#undef Sparse_Matrix
#endif